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Abstract 

The statistical properties of the E x B flux in different types of plasma turbulence 
simulations are investigated using probability density distribution functions (PDF). 
The physics included in the models ranges from two dimensional drift-wave turbu- 
lence to three dimensional MHD simulations. The PDFs of the flux surface averaged 
transport are in good agreement with an extreme value distribution (EVD). This is 
interpreted as a signature of localized vortical structures dominating the transport 
on a flux surface. Models accounting for a local input of energy and including the 
global evolution of background profiles show a more bursty behaviour in 2D as well 
as 3D simulations of interchange and drift type turbulence. Although the PDFs 
are not Gaussian, again no algebraic tails are observed. The tails are, in all cases 
considered, in excellent agreement with EVDs. 



1 Introduction 



In hot magnetised plasmas the main cross-field transport is often found to 
be far larger than expected for diffusive transport due to collisions. In plasma 
physics the transport is then called anomalous while in fluid dynamics " strange" 
or "anomalous" transport signifies the non-diffusive character of the transport. 
We investigate the statistical properties of the turbulent particle flux with the 
aim to examine the anomality of the transport in the fluid sense of meaning. 
An important feature stressed by many authors is the non-Gaussianty of the 
PDF of the point- wise measured flux [1,2,3,4]. While the non-Gaussianity of 
the flux PDF is to be expected from a quantity that itself is a folding of two sta- 
tistical variables, namely density and radial velocity, the observation of power 
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law tails in the flux PDF would be an indication of anomalous transport in 
the fluid sense of the wording, that is non-diffusive transport behaviour. Cor- 
respondingly this would indicate the presence of long range correlations in the 
turbulence [5]. Experimental results showing flnite size scaling and similarity 
of transport PDFs measured in different experimental devices have recently 
been reported in Refs. [6,7]. These measurements have shown that transport 
in the edge and scrape-off layer (SOL) region of fusion devices is highly in- 
termittent, self similar and does in its statistical features not vary much be- 
tween machines of different sizes and geometries. A characteristic feature of 
the transport seems to be the non-Gaussianity observed in its PDF and its 
intermittency. Intermittent transport would indicate that indeed the particle 
transport is anomalous not only in the sense that it is not due to collisions, but 
also in the sense that it is non-diffusive. A non-diffusive transport, however, 
has distinct scaling properties, which would pose much stricter conditions on 
the design of the vessel surrounding a fusion device, as this has to withstand 
the largest likely transport event. Additionally, a number of fundamental ques- 
tions in relation to transport PDFs are still unanswered, so it is as yet unclear 
what ingredients are needed to produce a certain transport PDF. Are trans- 
port PDFs dependent on the driving mechanism behind the turbulence? How 
does different physics enter into the transport PDF? Here we present initial 
answers to some of these questions by analysing direct numerical simulations 
of various types of plasma turbulence in different geometries. 



2 Local turbulence Models 

In this section we consider plasma turbulence in a classical sense, namely as 
being based on a scale separation between the equilibrium and fluctuations, 
f.x. density fluctuations are assumed small compared the the equilibrium den- 
sity n/noo << 1- The fluid equations for drift micro-turbulence result from 
standard ordering based upon the slowness of the dynamics compared to the 
ion gyro-frequency = eB/Mi and hence the smallness of the drift scale ps 
compared to the background pressure gradient scale length L±. These quan- 
tities and the sound speed Cg are defined by 

^'^W ''^S' ^^^S' ^^ = l^i°s^^^l"' W 

where subscripts e and i refer to electrons and ions respectively, and the tem- 
perature is given in units of energy. Normalization is in terms of scaled depen- 
dent variables (electrostatic potential e0/Te, electron density n/noo, parallel 
ion velocity u/cg, parallel electric current J/nooeCg; respectively). In addition 
the dependent quantities are scaled with the small parameter S = ps/L±, so 
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that wc calculate with quantities of order one. The quantity noo is a normaliz- 
ing density, while no{x) is the equilibrium plasma density having a finite gra- 
dient. In normalized units the radial profile of the density is dr logno(a;) = — 1. 
Vorticity is defined via a; = V^0. 



2.1 Two dimensional drift turbulence 



Assuming electrostatic motion, neglecting eletron inertia and replacing parallel 

derivatives by an effective parallel wavelength L,, we arrive at the well known 
Hasegawa-Wakatani equations [8] (HWE), describing two dimensional drift- 
wave turbulence in the absence of magnetic field line curvature: 

dtuj + {ip,uj} = --^ {n- (f) + fi^V'^uj, (2) 



dtn + {ip,n} + dyip ^ --^ {n - (f) + Hry'^n , (3) 



where the Poisson bracket if, o| = -J-if — irir is used to write the non- linear 

I.J T o J ax ay ay ax 

terms originating from advection with the E x B drift velocity. In this system 
the turbulence is driven by the resistive instabihty. We solve the HWE on a 
double periodic domain [9] and typically we use 256 x 256 modes on a square 
of side length L = 40. To produce self-consistent stationary turbulence we 
start the simulations from low amplitude (10~^) random density fiuctuations 
and initialise all other fields to zero. The fluctuations grow due to the resistive 
instability and saturate with amplitudes of order one. After saturation [10] we 
measure the local turbulent radial particle transport (x-direction) as given in 
a single point: 

^n,\0C = ' (4) 

and the flux surface averaged fiux T^pg: 

r„,Fs ^tI ^".loc dy ' (5) 

^ 



which in the two dimensional context is trivially obtained by averaging over 
the poloidal coordinate y. Density and potential fiuctuations have PDFs that 
are well described by Gaussians. Consequently the PDF of the E x B fiux, 
^nloc ^ depicted in Fig. 1, is close to the functional form expected from 
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the folding of Gaussian PDFs for the density fluctuations n and the radial 
E X velocity, Vr — —dy(j) [4]: 



^FG(r) 



exp 



(6) 



where c is the correlation between density and radial velocity fluctuations 
c = — , 2\i/jr/"2\i/2 — coscK^^^, Kq is the modified Bessel function of second 
kind, and anu is the phase angle between n and Vr- Furthermore we have 
a^^ = {v^y/'^{l — c^yl"^ , and correspondingly for Using Eq.(3) for an 
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Fig. 1. PDF of Jq^, compared with Fig. 2. PDF of F^pg, compared to an 
folded Gaussian EVD a log-normal distribution and a 

Gaussian 

approximate relationship between density and potential and the drift-wave 
dispersion relation to lowest orderin the long-wavelength limit u — j^^p- , we 
find: 



1 f LuU f 2 

^n,FS J ^^^^ '^'^^L J ^^y^^^'^ dy>0. 



(7) 



Thus the flux surfaced averaged transport is to lowest order positive definite, 
as confirmed by the numerical simulations. It is thus the PDF of the logarithm 
of the transport, InF^pg, that converges towards a Gaussian, if the central 
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limit theorem applies. As seen from Fig. 2 the transport PDF of F^^pg is 
indeed very well described by a log-normal distribution with in the present 

case zero offset z = 0: 




An alternative description of the statistical properties of the flux can be made 
in terms of the extreme value distribution (EVD). For a statistical variable - 
here the averaged flux - that is dominated by the minimum or maximum of a 
large number of random observations - here the local flux events - the EVD 
arises [11]: 



^EVD(r) = -exp 



exp 



exp 



9 , 



(9) 



The EVD has recently been used successfully to model and explain statistical 
properties of electron pressure fluctuations in a non- confinement plasma ex- 
periment [12]. As seen from Fig. 2 the EVD explains the observed transport 
PDF in this case as well as the lognormal distribution and is indistinguishable 
from the former. One should, however, note that the EVD decays faster than 
the lognormal one and, thus, is better suited to describe statistical variables 
that are limited not only from below, but also from above. Clearly the max- 
imum transport event in all real world systems is limited by the system size. 
We interpret the present situation in the sense that the flux surface averaged 
transport is dominated by extreme events, probably mediated by transport 
events linked to vortical structures. The flux PDF itself in this case is natu- 
rally skewed, but the radial particle transport is by no means 'strange' in the 
sense of not being diffusive on large time and space scales. No power law tail 
is observed in the data. That indeed the transport is diffusive is augmented 
by the fact, that the radial dispersion of ideal test-particles in 2D Hasegawa- 
Wakatani turbulence is found to be asymptotically diffusive as well [13,14]. 



2.2 Drift-Alfven and MHD turbulence 



Next we consider electromagnetic drift- Alfven turbulence with magnetic field 
curvature effects. The geometry is a three dimensional flux tube with local 
slab-like coordinates [a;, s] [15] and the model is deflned by the temporal 
evolution of the following four scalar flelds: electrostatic potential 0, density n 
given by electron density continuity equation, parallel current J, and parallel 
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ion velocity with the parallel magnetic potential A\\ — —V^J: 

— + Vb- Va; = X(n) + V||J + At^Via;, (10) 

+ Ve • V(no + n)^K{n-(f)) + Vy {J - u) + f^„Vln , (11) 

d 

— (Mil + fij) + fivE ■ VJ = V|| {no + n-(j))-CJ, (12) 
fdu \ 

el-^ + Vij-VMj =-V||(no + n) . (13) 



The advcctivc and parallel derivatives carry non-linearities entering through 
(f) and Ay: v^; ■ V = {(p, }, and Vy = |^ — }. Finally, the curvature oper- 

ator K = —ujb ismsdx + cossc?^^) originates from compressibility terms. The 
parameters in the equations reflect the competition between parallel and per- 
pendicular dynamics, reflected in the scale ratio e = {qR/Lj_)'^. The electron 
parallel dynamics is controlled by 

p^——e, i_i^—e, C = 0.51 H^i^H, 14 



where Tg is the electron coUision time and the factor 0.51 reflects the parallel 

resistivity [16]. The competition between these three parameters, representing 
magnetic induction, electron inertia, and resistive relaxation determines the 
response of J to the force imbalance in Eq. (12). The simulations are per- 
formed on a grid with 64 x 256 x 16 points and dimensions 64 x 256 x 27r in 
X, y, s. Some runs were repeated with double resolution to ensure convergence. 
Standard parameters for the runs were p, — 5, s — 1, ujb — 0.05 with the vis- 
cosities set to = fJ'n = 0.025. 

A transition in the dynamical behaviour from drift- Alfven to turbulence dom- 
inated by MHD modes is in this model triggered by increasing the ratio of 
P/jl. To check the nature of the observed turbulence we evaluate the phase 
shift a between density and potential fluctuations in dependence of ky using 
the relationship aky = log I {(f) kyUl^) . For better statistics the phase is aver- 
aged along the parallel coordinate s and time. Fig. 3 shows for a low value 
P = 0.5 and a high value /3 = 50 the change in principal behaviour of the 
turbulence with increasing (3. A shift towards larger phase angles is observed 
for low fcj^-modes indicating the transition from drift-Alfven to MHD type of 
turbulence. As the turbulence character changes to become more of the MHD 
type the phase relationship between density and potential is altered. This is 
also reflected in a visual inspection of the spatial transport structure (Fig. 4). 
It shows that the transport for high values of $ is carried by larger structures 
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Fig. 3. Probability of phase angle a be- 
tween density and potential fluctuations 
for /5 = 0.5 (top) and /3 = 50 (bottom) 
versus ky. 



Fig. 4. Spatial structure of the trans- 
port F^ Yoc 1^ ~ (left) and 
/? = 50 (right). 



than for low /3. The local flux F^ ^qq has the interesting structure that outward 
(white) and inward (dark) transport regions are closely related, this reflects 
the fact that the transport is carried in localized structures, having an outward 
flux at the one side and an inward flux at the other side. They are linked to 
corresponding eddies in the vorticity field. 

In both regimes the turbulence decreases the equilibrium density gradient by 
about 5 — 10% in the stationary turbulence, a further deviation from the 
background density gradient is prevented by a feed-back mechanism in two 
damping layers at the inner and outer radial boundary, which drive the flux- 
surface averaged density towards the initially specifled values. The PDFs of 
transport related quantities for /3 = 0.5 and /3 = 50 are shown in Figs. 5 and 
6. In the low beta case density and potential fluctuations have PDFs that are 
very well described by Gaussians, with minor deviations from Gaussianity in 
the tails. Density and potential have a rather similar PDF as they are well 
correlated and consequently the fltted folded Gaussian Ppc Eq. (6) describes 
rather well the point-wise flux PDF (c = 0.2) . An indication of a power-law 
tail in the point-wise flux PDF is not observed, and the point-wise flux PDF 
decays exponentially. The PDF of the magnetic flux surface averaged density 
flux F^ pg, shown in Fig. 5 c), is much closer to a shifted Gaussian, depicted for 
comparison in the flgure, but we observe a "fat-tail" towards larger transport 
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Fig. 5. Low (5 = 0.5 case, a) PDF's for 
density n and potential (j) fluctuations 
with fitted Gaussians (dashed) b) 
with fitted folded Gaussian c) F^ pg 
with a fitted Gaussian, log-normal and 
extreme value distribution (EVD). 
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Fig. 6. High /3 = 50 case, a) PDF's 
for density n and potential (j) fluctua- 
tions with fitted Gaussians (dashed) b) 
F^ loc '^ith fitted folded Gaussian and 
c) F^ pg with a fitted Gaussian and 
extreme value distribution (EVD). The 
log-normal distribution is not shown, as 
it is very close to the EVD. 



events. However, the flux surface averaged particle flux will, even though not 
positive deflnite, be at least limited from below, with excursions to a negative 
averaged flux being unlikely. Thus in Fig. 5 c) we also present the fltted log- 
normal distribution, with parameters a = 0.273, z = —9.3, and m = 11. This 
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fits the data well, but over-predicts the tail, an indication of the fact that the 
transport is limited from above as well, as we are considering a finite system. 
Indeed the EVD with parameters Hq = 2.60 and g = 0.34 fits the data very 
well, especially the tail of the PDF. For the MHD situation with ^ = 50 — 
even though the phase relation between density and potential differs and the 
turbulence has different character — we observe transport statistics similar to 
the drift-Alfvcn case. The corresponding PDFs for this case are depicted in 
Fig. 6. As the density and potential fiuctuations are less well coupled in that 
regime the PDFs of density and potential differ more. The potential fluctua- 
tions also deviate more from a Gaussian. Compared to the density fluctuations 
they now have a larger width. Note that the fluctuation level has increased 
significantly as well as the transport level when compared to the low beta 
case. Nevertheless the PDF of the point-wise measured transport is still well 
described by a folded Gaussian as seen in Fig. 6 b. The flux surface averaged 
transport is shown in comparison with a fltted EVD and a Gaussian. The log- 
normal distribution function is not plotted here, as it is very close to the EVD. 
The flux surface averaged transport behaviour has changed little as compared 
to the low beta case. It is still much closer to an EVD (or log-normal) than 
to a Gaussian one, however, we observe that to the low transport side the 
Gaussian distribution fits the transport PDF better. This effect is due to en- 
hanced levels of small noisy transport events occurring in the MHD type of 
simulation and the fiux PDF is more infiuenced by these small scale "random" 
transport events. Their distribution centers around zero, while the large scales 
reproduce the global transport PDF and carry the net fiux. The small scale 
transport events then make an appearance in the transport PDF for small av- 
eraged fiuxes, which explains the nearly Gaussian behaviour of the fiux PDF 
on the low transport side. 



3 Global Models 



After having considered local models we now drop the scale separation between 
turbulent fiuctuations and equilibrium — or in the cases where no equilibrium 
exists, as in the SOL, between fiuctuations and time averaged background. A 
motivation for this is that in the edge the fiuctuation level is of the same 
order of magnitude as the background and going out into the SOL region 
the fiuctuations in density can actually exceed the average density by a large 
factor. However, regarding fiuctuations and background on the same footing 
requires that the system is integrated to much longer times and thus higher 
demands in terms of computational power. We thus consider a two dimensional 
and a simple three dimensional model. 
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3. 1 Flute Model 



Concerning only two dimensional dynamics of the plasma is - if at all - only 
a good approximation in the SOL. We here consider a newly developed global 
interchange model for the full particle density ra, temperature T and vorticity 
with f 9, + 



djTi 

— + nK{4>) - K{n) = /x„Vln + - a„(n - 1) (15) 

O jrji J rji2 O 

2 at In 1 



((9t + i X V0 ■ V) a; - K{j)) = ii^V\uj - a^oj 



(17) 



Here the coupling between the equations is exclusively due to the curvature 
operator K = (d/dy where ( = 2ps/Ro with Rq the radius of curvature of the 
toroidal magnetic field. The coupling due to parallel currents giving rise to 
drift-wave dynamics is absent from these equations. The S'-terms on the right 

periodic 
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Fig. 7. Geometry showing shape of source and sink terms and the locations were 
long timeseries are taken. 

hand sides represent sources of particles and heat, while the a-terms operating 
in the SOL (Fig. 7) are sinks due to the parallel loss to end plates along open 
field lines. Boundary conditions correspond to free-slip, while in normalized 
units n = T= la.tx = Lri; and dn/dx = dT/dx = at a; = 0, see Fig. 7. 
For strong forcing the system is unstable to interchange modes causing signif- 
icant convective transport. A novel property of this model is the non-linear 
conservation of the energy integral E = J dx ^ (V0)^ -|- fnT , reveahng the 
conservative transfer of energy from thermal to kinetic form clue to magnetic 
field curvature. A realistic modeling of this process is mandatory for predictive 
global models. Nevertheless the solutions for strong forcing shows the charac- 
teristic "on-off" character of the turbulence due to self-sustained sheared flows 
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[17]. The time-avcragcd thermal gradients are strongest in the source regions 
and are "flapping" back and forth, leading to intermittent ejection of particles 
and heat far out into the SOL region. The transport PDF obtained from direct 





Fig. 8. PDF of pointwise and flux surface averaged transport inside LCFS and in 
the SOL 

numerical simulation of that model is very different from the ones obtained for 
the fluctuation based model. The flux surface averaged transport is no longer 
constant and instead varies with radial position as density is lost in the par- 
allel direction. We here take the transport data from a long timeseries since it 
is impossible to use spatial homogeneity in the radial direction to increase the 
quality of the statistics. It here makes sense to compare flux PDF's by normal- 
ising them to their root mean square. Fig. 8 shows the PDFs of the point-wise 
and the flux-surface averaged transport at radial positions corresponding to 
the edge and the SOL part of the simulation domain, (Pi and P4 in Fig. 7). 
A slight tendency to larger flux events in the SOL is observed. Averaging the 
transport poloidally has less dramatic effects than for the models considered 
in Section 2. The reason for this is the large poloidal correlation lengths of 
the order of the poloidal size of the simulation domain. Fig. 9 demonstrates 



1^ 0. 





Fig. 9. Tail of the flux surface averaged 
transport in the SOL (see Fig. 8) fitted 
with an EVD. 



Fig. 10. Flux surface and time aver- 
aged transport in the SOL fitted with 
an EVD and a power law. 



that even for this obviously highly intermittent and bursty transport the EVD 
{ho — 1.35 and g = 0.60 ) is an excellent fit for the tail of the PDF. However, 
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if the data are additionally averaged over times longer than the times between 
bursts, we observe a behaviour that can as well be approximated by a power 
law PDF as an EVD (Fig. 10). The amount of data being in the tail is small in 
that case, and the data-range is about 1.5 decades in F and down from six to 
only two decades in the ordinate (Fig. 9). One should note that whenever the 
variance of the data is large the body of a EVD will be similar to a straight 
line in a log- log plot, while its tail decays with increasing slope in the log- log 
plot. This shows the danger of fitting a PDF only for a very limited range in a 
doubly logarithmic plot. Our interpretation is thus, that to observe a distinct 
behaviour of the tail of the transport PDF one needs to be able to fit data over 
at least two decades and that fitting data on smaller regions is highly arbitrary 
and not suited to distinguish between different interpretations of the data. 



3.2 Electrostatic Drift Model 



Finally we consider a three dimensional global version of drift- wave equations. 
The model uses quasi-neutrality and electron density continuity, together with 
the parallel force balance equations for ions and electrons. Parallel convection 
is kept. No assumptions are made on the scale length of the background gradi- 
ent compared to the fluctuation scales, e.g. the constraint Vno ~ Vn ~ e <c; 1 
as underlying f.x. the HWE Eqs.: (2) and (3) is dropped. The resulting system 
is written in terms of the logarithm of density N — log(n): 



^ + VA^ • ^ V0 = fi^V^u; + V|| {U-V) + {U- l^) Vy AT (18) 

^ = - (vViiN + ViiV) + i^n (v^iV + (ViV)2) (19) 

^ = -V||0-l/e,^C/ (20) 

dV M 

- V|| (0 - iV) - (i/ei + l^e,n) V (21) 



dt m. 



e 



The energetic coupling is due to the parallel current dynamics that is by the 
difference in parallel electron velocity V and parallel ion velocity U. Now 
coupling due to curvature is absent from the equations. Again, the density is 
the full density evolving and is fed into the system on the inner side, while 
on the radially outer third of the system Bohm-sheath boundary conditions 
are implemented. The inner two thirds of the radius have periodic boundary 
conditions in the parallel direction. The y-direction is periodic and fixed values 
arc imposed on field quantities at the outer radial boundary, while the diffusive 
fluxes are set to zero at the inner radial boundary by prescribing a zero radial 
gradient there. The system is solved on a 128x64x30 grid. Due to the heavy 
computational load of that system it could not be integrated to times as high 
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as the 2D system, which influences the quahty of the statistics. The flux PDF 




Fig. 11. PDF of flux surface averaged transport in the edge (left) with fitted EVD 
on bulk and tail of distribution and in the SOL (right) with EVD fitted on the tail 
of the distribution. 

in this case shows distinct differences between edge region (closed field lines) 
and SOL with open field lines and losses to the limiter plates. For the edge 
region the poloidal correlation length is shorter than the box length and thus 
the PDF of the averaged transport is similar to that of the 2D- and 3D drift 
wave model. It fits well with an EVD in the tail and in the bulk (sec Fig. 11). 
The fiux PDF in the SOL part of the domain bears the characteristics of long 
poloidal correlations observed in that regime, due to the fiute modes - rather 
than drift waves - dominating the dynamics in the SOL. Thus the shape of the 
transport PDF, also after flux surface averaging, is very similar to that of the 
interchange model (see. Figs. 8 and 9), demonstrating also the applicability 
of 2D models in the SOL. However, the tail of the PDF again docs not show 
any sign of a power law behaviour and can over the whole range be fitted very 
well with an EVD. 



4 Conclusion 

Direct numerical simulations of physically distinct models of fiute and drift 
modes in two and three dimensional geometries have been presented and anal- 
ysed in terms of the transport PDF. For the fluctuation based models non- 
linear structures in the plasma dominate the transport, but give rise to a 
diffusive turbulent transport rather than a super-diffusion or otherwise anoma- 
lous (in the fluid sense) transport. Consequently a Gyro-Bohm like scaling of 
the transport with the magnetic held is found. Similar behaviour has recently 
been reported in investigations on test particle transport in plasma turbulence 
[14,18]. 

There is little difference in the results concerning different physical models. If 
the energy input into the system is changed from a homogeneous input due 
to a local, flxed pressure gradient to a localized plasma source and sink re- 
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gions are added the transport PDF becomes more bursty and longer poloidal 
correlations arise. Also in this situation the tails of all flux PDFs considered 
are extremely well fitted by extreme value distributions. Thus, the deviation 
of the point-wise as well as the flux surface averaged flux from Gaussianity is, 
at least in the asymptotic limit of large system size, not due to anomalous — 
in the sense of non-diffusive — behaviour. The presented observations should 
not be interpreted as indicators for the presence of self organized criticality. 
An interpretation in terms of occurrence of extreme transport events, which 
may be caused by locahzed eddies seems to be more appeahng. The statistics 
describing the transport PDF should therefor be taken from extreme value 
statistics, which in all cases flts the tails of the observed transport PDF very 
well. 
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